Imperfection effects for multiple applications of the quantum wavelet transform 
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We study analytically and numerically the effects of various imperfections in a quantum compu- 
tation of a simple dynamical model based on the Quantum Wavelet Transform (QWT) . The results 
for fidelity timescales, obtained for a large range of error amplitudes and number of qubits, imply 
that for static imperfections the threshold for fault-tolerant quantum computation is decreased by 
a few orders of magnitude compared to the case of random errors. 
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The mathematical theory of Wavelet Transforms (WT) 
finds nowadays an enormous success in various fields 
of science and technology, including treatment of large 
databases, data and image compression, signal process- 
ing, telecommunications and many other applications 
mm Wavelets are obtained by translations and dila- 
tions of an original function and they allow to obtain 
high resolutions of microscopic details, both in frequency 
and space. The discrete WT can be implemented with 
high computational efficiency and provide a powerful tool 
for treatment of digital data. It is well accepted that the 
Fourier transform and WT are the main instruments for 
data treatment, and it has been shown that in many ap- 
plications the performance of WT is much higher com- 
pared to the Fourier analysis. The permanent growth of 
computer capacity has significantly increased the impor- 
tance of the above transformations in numerical applica- 
tions. 

The recent development of quantum information pro- 
cessing has shown that computers based on laws of quan- 
tum mechanics can perform certain tasks exponentially 
faster than any known classical computational algorithms 
(see e.g. Q). The most known example is the integer 
factorization algorithm proposed by Shor 4]. An essen- 
tial element of this algorithm is the Quantum Fourier 
Transform (QFT) which can be performed for a vector 
of size N = 2 ,lq in 0{n£) quantum gates, in contrast 
to 0(2 nq n q ) classical operations 0. EJ Here n q can be 
viewed as the number of qubits (two-level quantum sys- 
tems) of which a quantum computer is built. Apart from 
Shor's algorithm, the QFT finds a number of various ap- 
plications in quantum computation, including the simula- 
tion of quantum chaos models showing rich and complex 
dynamics 0, 0|- The sensitivity of the QFT to im- 
perfections was tested in numerical simulations and the 
time-scales for reliable computation of the algorithm were 
established [S 011 El 

A few years after the discovery of the QFT algorithm, 
it has been shown that certain WT can also be imple- 
mented on a quantum computer in a polynomial number 
of quantum gates [HI HH E2] In fact; explicit quan- 
tum circuits were developed for the most popular discrete 



WT, namely the 4-coefficient Daubechies WT (D^) and 
the Haar WT, both for pyramidal and packet algorithms 
[lCt 111! [12J . As it happens in classical signal analysis, it 
is natural to expect that QWT will find important future 
applications for the treatment of quantum databases and 
quantum data compression. Therefore, it is important to 
investigate the stability and the accuracy of QWT in re- 
spect to imperfections. This is especially important since 
the functions of the wavelet basis have singularities in the 
derivatives (in contrast to analyticity of Fourier waves) 
that may enhance the effects of perturbations. 

To this aim we introduce a simple model with rich 
nontrivial dynamics which is essentially based on multi- 
ple applications of the WT. Its quantum evolution can 
be efficiently simulated on a quantum computer, and it 
is described by the unitary map for the wave function ip: 



(1) 



Here the bar marks the value of the wave function af- 
ter one map iteration, W is the WT operator, 
and the unitary diagonal operators Ut = e~ lTn / 2 and 



-ik(x--rr) 2 /2 



represent quantum phase rotation 



in computational and wavelet basis, respectively. The 
evolution takes place in the Hilbert space of N — 2 rlq 
states, with —N/2 < n < N/2 and x — 2nj/N where 
j = 0, ... N — 1 is the index in the wavelet basis and T, 
k are dimensionless parameters. In the case when W is 
replaced by the Fourier transform, one obtains the quan- 
tum sawtooth map previously analyzed in Ref.Q. Thus 
the model Q can be considered as a 'kicked wavelet ro- 
tor', where k is the kick strength in the wavelet basis. 
We numerically tested that the dynamical properties are 
not very sensitive to the value of T and here we present 
data for a typical value T = 1.4. 

The global properties of the evolution operator (JJJ are 
shown in Figs lll2l for different values of k (see also Ap- 
pendix Al). The density plot of transition matrix el- 
ements U n ^ n i in the computational basis is represented 
in Fig^ By increasing k a larger and larger number of 
states is coupled by the dynamics, and the complex self- 
similar structure of the transitions generated by the WT 
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FIG. 1: (Color on line) Density plot of matrix elements 
\U n ,n'\ 2 for the model Q in the computational basis, for 
N = 2 12 . Top: k = 100 (left), k = 1000 (right); bottom is for 
k = 1000: a doubled resolution of left upper quarter (left), 
perturbed operator with static errors e = 5 x 10- 4 ,fi = 0. 
Color marks the density from blue (zero) to red (maximal 
value) . 



becomes evident. On the average, the off-diagonal matrix 
elements decay with the power law | | 2 ~ l/\n— n'\ a . 
Asymptotically for \n — n'\ 3> 5k we obtain the exponent 
a — 4 (FigEJl. For large values of k the intermediate 
scaling law is described by the exponent a = 2, in the 
range 1 < \n — n'\ <C 5k. This decay law for the ma- 
trix elements can be considered as a long range coupling 
between states. We note that similar power law reg imes 
have been analyzed in random matrix models |l3l 1 1 - ll ] . 
Our numerical analysis shows that there are two regimes 
for the level spacing statistics P(s) 15} in the limit of 
large TV. E.g. for N — 2 12 the distribution P(s) is given 
by the Poisson law for k < 5, while for 5 < k < 10000 it 
shows level repulsion and a poissonian decay for large s 
(see Appendix A2). We attribute the rapid appearance 
of level repulsion to the slow power law decay of matrix 
elements |l4| . 

To implement the evolution (J on a quantum com- 
puter, we developed an algorithm based on the QWT 
for the Daubeschies D 1 - 4 ^ wavelets. The algorithm con- 
sists of four steps: i) the multiplication by Ut, performed 
in 0(n q ) controlled-phase shift gates as described in [2|; 
ii) the application of W operator, realized by the QWT 
following the circuit described in Fig. 10 of [llj (see Ap- 
pendix A3); iii) the operator Uk, implemented in a similar 
way as for the step i); iv)the inverse WT W\ obtained 
by reversing the gates of the step ii). The heaviest parts 
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FIG. 2: Dependence of averaged matrix elements (\U n ^ n i j 2 } 
on \n — n'\ (the average is taken along the diagonal). Data 
are shown for N = 2 15 and k — 1 (full black curve), k = 10 
(dash-dotted curve), k — 100, (dotted curve) and k = 1000 
(dashed curve). The two straight lines are l/|n — n'| 2 and 
l/\n . — n'| 4 . The inset shows the data in semi- log scale. 



of the algorithm are the steps ii), iv), since the QWT al- 
gorithm requires multi-controlled operations. To imple- 
ment them we used the recipe given in [lfij | which allows 
to realize a n-controlled gate by 0(n) elementary gates 
(Toffoli and 1- and 2-qubit gates). To this end an ancilla 
qubit is needed, so that we used n q + 1 qubits to simu- 
late numerically the dynamics of model Q with N = 2"« 

states. The implementation of the wavelet kernel D$ re- 
quires 0(n) multi-controlled gates (n = 2, . . . ,n q ), and 
since the QWT is composed of 0(n q ) kernel applications 
this leads the total number of elementary gates to scale as 
0{n\) 0EJO3 (see Appendix A3). To study the algo- 
rithm accuracy we consider two models of imperfections. 
In the model of random noisy gates we replace all ideal 
gates by imperfect ones, which are obtained by random 
unitary rotations by a small angle 77, — e/2 < r/ < e/2, 
around the ideal rotation angle (as in 01)- I n * ne model 
of static imperfections (see |j1LL£|) all gates are perfect 
but between gates ip accumulates a phase factor e 1 ^ with 
4> = ^liVlVi+Mvfvi+i)- Here r/i, fi t vary randomly with 
I = 0,...,n q , rji represents static one-qubit energy shifts, 
—e/2 < rji < e/2, and /i; represents static inter-qubit 
couplings on a circular chain, — /i/2 < /i; < /i/2. 

The numerical simulations of the ideal quantum algo- 
rithm for the map Q show that the wave function is 
essentially localized on a few states of the computational 
basis. This localization is clearly seen from the Inverse 
Participation Ratio (IPR) £ = 1/Xml^™l 4 wn i cn i s a 
standard quantity to characterize localization in meso- 
scopic systems [Taj . It directly provides the number of 
sites on which the probability is concentrated. Surpris- 
ingly the localization is present not only for moderate 
k ~ 1, but also when the kick strength is very large 
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FIG. 3: (Color on line) Dependence of IPR £ on the number 
of iterations t, for n q = 12, T = 1.4, fc = 1 (top) and k = 1000 
(bottom). Initially the probability is concentrated at n = 0. 
The black curves show the quantum computation with ideal 
gates; the green (light gray) curves show the case with static 
errors at e = 10~ 4 ,^i = and red (gray) curves correspond 
to the case with noisy gates at e = 5 x 10~ 4 . The data are 
averaged over time interval At = 50. 



k ~ 1000 (see Fig. [3 and Appendix A4). Indeed in both 
cases £ fluctuates near a constant value £o N, even for 
a very large number of iterations. We attribute this local- 
ization to the structure of the operator (JIJ: it is banded 
for moderate k and sparse for large k (see Fig.Q). For 
k ~ 1 the probability shows an algebraic localization 
\ip n \ 2 oc 1/n 4 (FigQJ. Such an exponent fully agrees 
with the scaling law of Fig|21 For k > 100, the prob- 
ability is spread over the whole basis (data not shown), 
but only a moderate number of narrow peaks contributes 
to the IPR value (see Fig0 ■ This behaviour is consistent 
with the fact that the P(s) never reaches a Wigner-Dyson 
regime (see discussion above) . On the contrary, the spec- 
tral properties of the sawtooth map 0, E| are described 
by the random matrix theory for k ~ 1000, T ~ 1 and 
N = 2 12 . 

The effect of imperfections in the quantum gates is 
shown in Figs. 01El( see a ls° Appendix A4). The results 
clearly show that the localization is destroyed by noisy 
gates imperfections which lead to an approximately lin- 
ear growth of £ with t. For static imperfections £ shows 
modified bounded oscillations. The probability distri- 
bution in Fig0] shows the appearance of a plateau with 
pronounced peaks located approximately at n = JV/2™, 
m = 1,2,3.... We attribute the appearance of these 
peaks to the pyramidal structure of the algorithm, which 
in the presence of imperfections produces stronger er- 
rors at the above values of n. For static imperfections 
the plateau level remains bounded in time t while for 
noisy gates it increases with t and for very large t the 
probability becomes homogeneously distributed over the 
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FIG. 4: Probability distribution \ipn\ 2 in the computational 
basis for the parameters of Fig0 (top) at k = l,t = 10 4 : 
full curve is the quantum computation with ideal gates, gray 
curve shows data for noisy gates with e = 5 x 10~ 4 . The 
straight line displays the scaling law 1/n 4 . The inset shows 
the same data in semilogarithmic scale. 

computational basis. 

The qualitative difference between two types of imper- 
fections becomes clear from the analysis of the behaviour 
of the fidelity, defined as f(t) = \(ip e (t)\ip(t))\ 2 . Here tp(t) 
is the wave function obtained with ideal gates, while ip e (t) 
is the result of the quantum computation with imperfec- 
tions of amplitude e. We determine the time scale tf 
for accurate computation by fixing a threshold for the 
fidelity as f(tf) — 0.9. In this way it is possible to find 
the dependence of tf on the system parameters. Our nu- 
merical data are presented in FigOO They show that for 
noisy gates tf is described by the relation 

t f = C/(e 2 n g ), N g = C/e 2 (2) 

where n g is the number of gates per map iteration, N g — 
n g tf is the total number of gates and C ~ 5 is a numerical 
constant. The physical origin for this scaling is related 
to the fact that after each gate an amount of probability 
of the order of e 2 is transferred from the ideal state to all 
other states. This leads to an exponential decay of the 
fidelity f(t) w exp(— Ae 2 n g t), where A is a constant (see 
FigEt). This gives the scaling which was also found 
in other algorithms with noisy gates 0, El El • 

For the model with static imperfections the scaling is 

t f = D/(en g n q 1 / 2 ), N g = D /{en, 1 / 2 ) (3) 

where D is a numerical constant (_D«4.5, at/i = and 
D w 2.1 at fi = e). This timescale is significantly smaller 
than the one for noisy gates. Physically, this happens due 
to the coherent action of static imperfections, which lead 
to effective Rabi oscillations proportional to cos (en g t) for 
each qubit. For n q qubits this gives f(t) cx [cos (en g t)] nq 
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FIG. 5: Panel a) shows the fidelity decay in time at k = 1, 
n q = 12 for static imperfections (e = 10 _4 ,^i = 0, dashed 
curve) and noisy gates (e = 5 x 10~ 4 , full curve). Panel b) 
shows the dependence of time scale £/ on the imperfection 
strength e for n q = 8 (n g — 5237) for noisy gates (diamonds) 
and static imperfections (triangles at n = 0; circles at fi — e, 
for clarity data are shifted in e-axis by factor 10 to the left). 
Panel c) gives the dependence of the total number of gates N g 
on e for n q = 6, 8, 10. For noisy gates (diamonds) e = e and for 
static imperfections (triangles) e = e^/n^. Open(full) symbols 
are data for k — l(k — 1000). The full and dashed lines in 
panels b)/c) show the relations and ©, respectively. 

and for small e we obtain a Gaussian drop of the fidelity 
f(t) ~ exp (— n q (en g t) 2 ), in agreement with our numer- 
ical results (see Fig|^). This leads to the scaling |JjJJ, 
which is confirmed by the data in Fig|S] The effects of 
static imperfections are dominant for all range of imper- 
fection strengths studied. We note that the scaling laws 
@ , ® are rather general and do not depend on the kick 
strength k. We note that similar scalings were discussed 
and numerically demonstrated in other quantum algo- 
rithms with noisy gates 0, Il7| and static imperfections 
(see also HJj)- This shows that such scaling laws are 
generic and are not sensitive to the singularities in the 
derivatives of the wavelets. The universality of the above 
relations iffll. iffll is also confirmed by the fact that the 



structure of the QWT is rather different from the QFT 
algorithm, e.g. the number of elementary quantum gates 
scales as 0(n^) for the QWT, in contrast to 0(n 2 ) for 
the QFT. These relations determine the total number of 
gates N g — tfUg during which the quantum computa- 
tion is reliable. Similar scalings for N g should also be 
valid for other quantum algorithms, e.g. Grover's and 
Shor's algorithms. We discuss also other types of errors 
in Appendix A5. 

The above relations J2J, © are important for the 
quantum error correction codes and the fault-tolerant 
quantum computation threshold (see 0, and Refs. 
therein) . Indeed the accuracy border for large scale quan- 
tum computation is obtained in the assumptions of ran- 
dom noisy errors and gives a threshold e < e r ~ 10~ 2 . 
This approach intrinsically uses the fact that for noisy 
gates the fidelity remains close to one for a number of 
gates N g — C/ef (see (0)). In the case of static im- 
perfections it is natural to assume that this number of 
gates should remain approximately the same to allow 
large scale computation on a quantum computer with 
n q qubits. Therefore, for static imperfections the Eqs. 
(0), 10 give the accuracy border e s : 

e s « De 2 r /(Cn^ 2 ) (4) 

This important relation gives a significant decrease of the 
threshold for the case of static imperfections [z| • For the 
parameters of our model at n q — 10 we obtain that for 
the noisy error rate p r — e 2 . « 10~ 4 the rate induced by 
static imperfections should be less than p s = e 2 s w 10 -9 . 
This result shows that new strategies of quantum error 
correction codes should be developed to significantly sup- 
press phase shifts induced by static imperfections. The 
spin echo techniques used in NMR [2| may play here an 
important role. 

This work was supported in part by the EC contracts 
RTN QTRANS and IST-FET EDIQIP and the NSA and 
ARDA under ARO contract No. DAAD19-01-1-0553. 
We thank CalMiP in Toulouse and IDRIS at Orsay for 
access to their supercomputers. 
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APPENDIX 
Al 



Here we show examples of density plot for the matrix elements |[/ Tli „<| 2 for the model (QJ for k — 1, 10, 100, 1000 
(Fig©- 




FIG. 6: Density plot of matrix elements |t/„ in '| 2 for the model Q in the computational basis, for TV = 2 12 , k = 1 (top left), 
k = 10 (top right), k — 100 (bottom left) and k = 1000 (bottom right). Color marks the probability density, from blue to red 
(maximal value). 
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A2 

The spectral analysis of the model is obtained by a numerical diagonalization of the evolution operator U. 
Due to the unitarity of U, all the eigenvalues A are on the unitary circle, A = e luJ , where u> are the quasi-energies 
included in the interval [0, 2n). The typical examples for the level spacing statistics P(s) for uj are shown in FigH 
It is remarkable that even for large kick strengths {e.g. k — 1000) the Wigner-Dyson statistics of the random matrix 
theory is not achieved. We note that for such values of k in the sawtooth map all the eigenstates are delocalized and 
P(s) is given by the Wigner-Dyson distribution 

0E3. 
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FIG. 7: Level spacing statistics P(s) for the quasi-energies of model |T), for n q — 12 and different values of the parameter k. 
A transition from the Poisson distribution P(s) = e~ s at small k values to a distribution which shows the level repulsion for 
small s is observed by increasing k. Data are shown for k = 0.1 (squares), k = 1 (circles), k = 10 (diamonds) and k = 1000 
(x's). The inset displays the data for k = 1 and k — 10 in a semilogarithmic scale. The full line is the Poisson distribution. 
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A3 

Our quantum circuit is based on the scheme described in We implemented the Pyramidal Algorithm (PYA) 
for the Z?( 4 ) wavelet transform. It is based on repeated applications of the operator D^} (the wavelet kernel) and 
the permutation operator n 2 ™. Here the kernel D$ is the Daubeschies matrix of size 2™ x 2™. The operator 
II 2 n realizes the shuffling step on vectors {u?}j=i,2 n of size 2™. The action of H 2 n can be regarded as a (classical) 
permutation of the index for the vector {i>j}j=i,2 n > The binary representation of index j, (ao, a,\, . . . , a„_i), is mapped 
into (a„_i, ao, a\, ■ ■ ■ , a n ^ 2 ), where ao is the most significant bit. The classical operator can be written as 

= {D<£ ] © 7 2 n 9 -4)(n 8 © / a »,_ 8 ) . . . {Df} © / a «,_ 2 *)(n a * +1 © J 2 » 9 _2*+0 . . . n 2 »,D^ (5) 

where ijvf is the identity matrix of size M x M and © is the direct sum of operators (see FigJHJ) . 

In a quantum computation, the action of Il2« on the element \j) = \ao,a±, . . . ,a„_i) of the computational basis 
is Il2» |eto, (2i, . . . , a n —2, a„_i) = |a„_x, a n -2, ■ ■ ■ ,a>u a o) an d it can be implemented via n — 1 quantum swaps, each of 
them built by 3 control-not gates. The direct sums 112" ®I 2 nq- 2 n and Z) 2 t ©/2™9-2™ correspond to multi-controlled 
operators with n q — n controlling qubits. Both II2™ and D$ can be implemented by a polynomial sequence of 
elementary gates (Toffoli, Control-Not, and one-qubit rotations), therefore the above multi-controlled operators are 
replaced by the product of multi-controlled elementary gates. Following the procedure proposed in [16j, these multi- 
controlled gates were implemented through elementary gates with the help of an ancilla qubit. The computational 
cost of a ^-controlled gate is linear in the number of controlling qubits I. 

The wavelet kernel D$ is decomposed into elementary gates following the factorization proposed in with slight 
modifications. The kernel can be written as 

£> 2 n = (ho-l <8> Cl)P 2 »(iV ® I 2 n-l)(N © I 2 n-2 © I 2 n-x) ... (N © I 2 © I 2 ~-i){N © I 2 n _ 2 )P 2 „ (I 2 n-1 © C ) (6) 

where P 2 ™ is the full permutation matrix which action on the binary representation of vector indexes is 
P 2 n(a,Q,ai, . . . ,a n - 2 ,a n —i) = (a n _i, a n _ 2 , . . . , ai, ao) and N is the not gate. Here C\, Cq are 2x2 rotation ma- 
trices, which can be expressed via the Daubechies coefficients cq, c\, c 2 , C3 by defining 

*-S(?4) 

Q _ 1 ^ _ f sin O cos O \ C — ^ ^ _ / sin6>i cos6>i 

° ~ VdetCo ° V cos0 o -sin0 o J 1 ~ ^dotCx" 1 \cos0i ~ sin ^ 

where $o = § an d 9\ = j^tt. The operator P 2 ™ is implemented by 0{n) swap gates. We note a slight modification 
in the equation J7J, comparing to The quantum circuit corresponding to the wavelet kernel © is shown in 

Fig|5] Fig ^| clarifies the notations used in Figs ECU The total number of elementary gates needed to implement the 
kernel circuit scales as 0(n 2 ), thus leading to a 0(n q ) total complexity for the QWT. For our circuit the number of 
elementary gates was n g = 1509, 2974, 5237, 8470, 12821, 18462, 25541 for n q = 6, 7, 8, 9, 10, 11, 12. It is assumed that 
the elementary gates act between any two qubits. 
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FIG. 8: Quantum circuit for the wavelet Transform @. 
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FIG. 10: Representation of the (A 2 n-i © I 2 n-i) operator, @ is the Not Operation. 
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A4 

Here we show the probability distribution in the computational basis for a large value of kick strength (k = 1000) at 
two different moments of time (t = 1000 and t = 10000) (Figs lllJTSjl . We remark two main features: the distributions 
have pronounced peaked structure and the peaks are spread all over the computational basis. For t = 1000 the effects 
of noisy errors are weak so that exact and noisy distributions are close (top vs. middle) while the distribution for 
static imperfections is already strongly modified (top vs. bottom). 
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FIG. 11: Probability distribution in the computational basis for n q = 12, k = 1000 and t = 1000: quantum computation with 
exact gates (top), with noisy gates at e = 5 x 10~ 4 (middle) and static imperfections at e = 10~ 4 ,fi = (bottom). 
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FIG. 12: Same as Fig[TT]for t = 10000. 
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A5 

We have also considered another model of static imperfections. It is obtained from the model of noisy gates by 
repeating the same sequence of errors for each application of the evolution operator U in Q). As in ^3] each noisy 
gate transformation is obtained by diagonalization of nondiagonal part and then by multiplication of each eigenvalue 
by a random phase exp(z?y) with — e/2 < r\ < e/2. This pseudo-static imperfections model is intermediate between 
the two cases considered in the text. The behaviour of fidelity f(t) is similar to the one shown in Fig|S^ (see Figll3l 
(top)). For large and moderate e the total number of gates N g is not very large compared to n g and correlations 
between different map iterations can be neglected. Then the gates look like quasi-random and the scaling is given by 
the relation (J2J with C » 5. However, in the limit of small e the coherent rotations become dominant and the data 
give the scaling N g ocl/e (see Fig llHJI . This confirms the generic scaling typical of static imperfections. 
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FIG. 13: Comparison between noisy gates model and pseudo-static imperfections model. Top panel shows the behaviour of 
the fidelity f{t) for n q = 12, k = 1, T = 1.4 and e = 5 x 10~ 4 (full curve for noisy gates model) and e = 10 4 (dashed curve 
for pseudo-static imperfections). Bottom panel: scaling for the total number of gates N g as a function of e for noisy gates 
(diamonds) and pseudo-static imperfections (triangles). Open (full) symbols correspond to k = 1 (k = 1000). The full and 
dashed straight lines show the dependences N g oc 1/e 2 and N g ocl/e for noisy gates and psedo-static imperfections respectively. 
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